A Theory for the Membrane Potential of Living Cells 
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We give an explicit formula for the membrane potential of cells in terms of the intracellular 
and extracellular ionic concentrations, and derive equations for the ionic currents that flow through 
channels, exchangers and electrogenic pumps. We demonstrate that the work done by the pumps 
equals the change in potential energy of the cell, plus the energy lost in downhill ionic fluxes through 
the channels and exchangers. The theory is illustrated in a simple model of spontaneously active 
cells in the cardiac pacemaker. The model predicts the experimentally observed intracellular ionic 
concentration of potassium, calcium, and sodium. Likewise the shapes of the simulated action 
potential and five membrane currents are in good agreement with experiments. We do not see any 
fZ—i . drift in the values of the concentrations in a long time simulation, and we obtain the same asymptotic 

' values when starting from the full equilibrium situation with equal intracellular and extracellular 

ionic concentrations. 

^! I. INTRODUCTION 

o 

The purpose of the work we present here is to obtain a model for the membrane potential of a single cell which is 
reasonably realistic, and yet so simple that it can be used in practice to simulate numerically single cells or several 
coupled cells. For simplicity, experimentally observed currents (Boyett et al. 1993) are omitted if they seem too small 
r/} \ to have a significant influence on the intracellular ion concentration, or at least too small to change the dynamics of 
the cell. On the other hand, we try to make the theory realistic by using equations that are compatible with, or can 
be derived from, basic physical principles. 

It is a basic assumption of most models (Wilders 1993) for the electrical activity of cells that only the motion 
of positive ions, and specifically those of potassium, calcium and sodium, influence the membrane potential. This 
assumption is usually expressed as a differential equation for the time dependence of the potential. We observe that 
this differential equation can be integrated exactly, and argue that the integration constant is given by the requirement 
I ■ that the potential is zero when the ion concentrations on both sides of the membrane are equal, as the density of 
negative charge happens to be the same on both sides. Then it follows that the potential is directly proportional to 
t-H . the excess number of positive ions inside the cell, a formula which is nothing but the one for an electric capacitance 
that follows from Gauss's law in electrostatics. 

We derive equations for ionic currents flowing through channels, exchangers and electrogenic pumps. These are 
based on the Boltzmann distribution law (Boltzmann 1868), which states that a particle in thermal equilibrium spends 
less time in states of higher energy than in states of lower energy, the Markov assumption (Markov 1906) which says 
• i-H | that the transition probabilities of a stochastic system (of Markov type) is only dependent on its present state, and 
the principle of detailed balance (Onsager 1931) which says that the microscopic laws of physics are invariant with 
respect to the reversal of time. Our equations were inspired by Ehrenstein and Lecar's model of channel gating (1977), 
Q_i Nonner and Eisenberg's model for channel current (1998), Mullins' model of the Na + ,Ca 2+ exchanger (1977), and 
Chapman's model of the Na + , K + pump (1978). In particular the book of Lorin John Mullins (1981) "Ion Transport 
in Heart" has been a major source of inspiration for us. 

The theory is illustrated with a simple model of spontaneously active cells in the rabbit sinoatrial node. The 
observable parameters in the model are based on the experiments of Shibasaki (1987), Hagiwara et al. (1988), 
Muramatsu et al. (1996) and Sakai et al. (1996). The non-observable parameters in the model are determined 
numerically, in the same way as in an earlier study (Endresen 1997a), by comparing the action potentials generated 
by the model with the shape of the action potentials recorded by Baruscotti et al. (1996). 

By using an algebraic equation for the potential in place of the standard differential equation, as mentioned above, 
we obtain a model which is stable against a slow drift of the intracellular ion concentrations, sometimes seen in other 
models. Furthermore, by fixing the integration constant for the voltage we obtain from the model a prediction of the 
steady state ion concentrations in the cell. It is even possible to predict these steady state concentrations by starting 
with an initial state having equal concentrations inside and outside the cell, and integrating the equations of motion 
over a long time interval. 
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From the equations of motion we obtain an equation that explicitly demonstrates the energy balance in the process 
of moving ions in and out of the cell. The energy to make the cell function comes from breakdown of ATP that runs 
the Na + , K + pump. Part of this free (or useful) energy is dissipated while the rest enters the cell. In the cell some of 
this energy is used to create a potential energy that depends upon the ionic concentrations in the cell, while the rest 
is dissipated by the currents in the ionic channels and the Na + ,Ca 2+ exchanger. This potential energy function is 
thus such that the work associated with ionic currents balances exactly the change in potential energy. In a numerical 
integration of the differential equations one may compute separately the work and potential energy, comparing the 
two in order to check (and maybe control) the accuracy of the numerical integration. In our long time integration we 
observe indeed a balance between work and change in cell membrane potential energy. 



II. DERIVATION OF THE EQUATIONS 

A. The Nernst Equilibrium Potential 

There are two basic principles behind the average motion of particles. The first is diffusion, which is general; the 
second applies only to charged particles such as ions in solutions. Simple diffusion is described by the empirical law 
of Fick (1855), 

cf=-ukTV[S], (1) 

where <fi is the ionic flux, [S] the concentration of ions and u the ratio of the velocity to the force acting on a particle, 
known as the mobility. T is the absolute temperature and k is Boltzmann's constant. The empirical law of Ohm 
(1827) describes the net motion of charged particles in an electric field, 

$ = -zeu[S]\7U , (2) 

where z is the valence, e the elementary charge and U the electrical potential. Since we assume that the temperature 
is constant, we can neglect the thermal flux given by Fourier's empirical law. The fact that the mobility in Fick's law 
must be identical to the mobility in Ohm's law was first noticed by Einstein (1905). If we combine Eqs. ([!]) and <^), 
the total flux of ions due to diffusion and electric forces is 



The equilibrium potential for which the flux is zero, is 



[S] exp 



zeU 



(3) 



v s = U i -U e = ^ lnfgh . (4) 



[S]i 



It can be found by setting <ft — in Eq. (Q) and integrating from the extracellular (e) to the intracellular (i) side 
of the membrane. Here Ui, U c , [S]i and [S] c are the intracellular and extracellular potentials and concentrations. 
This equation, first stated by Nernst (1888) is based only on the empirical laws of Ohm and Fick and the relation of 
Einstein. 

The same formula can be derived in a more general way using the Boltzmann factor (Boltzmann 1868). The relative 
probability at equilibrium that an ion is at the intracellular or extracellular side of a cell membrane is 

Pi [S], / ze(U,-U c ) \ 

where ze{U\ — U e ) is the energy difference between the two positions of the ion. Solving (§) for U - U c gives (g). 
With the definition 

kT RT 

vt= — = , 6 
e t 

the equilibrium potentials for the predominant cellular cations are then 
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(7) 
(8) 
(9) 



B. Ionic Channels 

i. Ionic Channel Gating 

Imagine that ionic channels are either completely open or completely closed and randomly fluctuate between these 
states in a simple Markov process (Markov 1906), described by the first order kinetics (Ehrenstein and Lecar 1977) 

where the rate constants a and are functions of transmembrane voltage and control the transitions between the 
closed (G) and the open (O) states of the gate. The rate for a closed channel to open is a, and (3 is the rate for an 
open channel to close. Let x denote the average fraction of channels that are open, or, equivalently, the probability 
that a given channel will be open. We may say that the ionic flux through an ensemble of channels is regulated by a 
sliding door whose position is x. This yields: 



df 
where 



a (i _ x ) _ p x = ^_J^ f (n) 



a 



a + P 
1 

a + P 



(12) 
(13) 



Here Xqq denotes the steady state fraction of open channels and r the relaxation time. Let us assume that the energy 
difference between the open and closed positions is given by 

AG = G opon - G c i osc d = q(v K - v) , (14) 

where q is a gating charge, usually q w ±4e, such that qv represents the change in electrical potential energy due 
to the redistribution of charge during the transition, and where the term qv x represents the difference in mechanical 
conformational energy between the two states. At equilibrium, dx/dt = 0, and the ratio of the probabilities for a 
single channel to be in the open state or the closed state is 

= a (lg) 

This relation is known as the principle of detailed balance (Onsager, 1931). The same ratio is given by the Boltzmann 
distribution (Boltzmann 1868), 

( AG\ , . 

ex P ( ~Hm ) • (16) 



Thus, from Eqs. (||), (jTJ), and (UJ), with q = +4e 



exp 



Ae(v x - v) 
kT 



(17) 



The simplest possible choice for a and (3 is 
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a = Aexp 
= A exp 



2e(t> x - ^) 

fcT 
2e(i; x - t>) 

fcT 



(18) 
(19) 



where A is a constant. Assuming A to be constant gives the maximum relaxation time at the voltage where Xoo = 1/2. 
The relaxation time as a function of v is then 



1 

a + /3 



2A cosh 



2e(v x - v) 
kT 



(20) 



2. Ion Channel Current 



Here we want to obtain the current i through a one-dimensional ionic channel of length d. We will allow the cross 
sectional area A to vary with position, i.e., we take A = A(x). By definition, x = —d/2 is the inside and x = d/2 
the outside of the membrane. Let <fi = 4>(x) be the x-component of the flux <fr, the other components are negligible as 
long as the variation of A with x is smooth and slow. This is the analogue of water flow in a pipe of varying cross 
section. By stationary flow, the current i must be the same through all cross sections, i.e. independent of x. Thus 
the flux (j) is inversely proportional to the area A, by the relation 



i = ze(f>A = const. 



(21) 



We insert 4> from this equation in the x component of Eq. (|^), and multiply the resulting equation by 
exp (ze(U — Uo)/kT), introducing a constant voltage Uq chosen such that 



U 



U 



Then we find that 



/ ze(U - Uq) 
exp ' 



.4 



kT 



-zeukT 



U 



U 



dx 



[S] exp 



ze{U - U ) 
kT 



(22) 



(23) 



Here U, [S] and A are functions of x, while all other quantities are constant. (Note however that the mobility u may 
be reduced in a very narrow channel; one possible way to take into account such an x dependence of u is to replace 
the area A by an effective area A c $ which is smaller than A). Integrating from the inside x — —d/2 to the outside 
x = d/2 we obtain 



ikT 



where 



[sl . exp ,_-)_ [S Mli) 



d/2 



.4 



kT 



(24) 



(25) 



The concentrations are [S]; on the inside and [S] on the outside. If we extract a factor y^S]7[S]i, and write the ratio 
of the concentrations in terms of the Nernst potential defined in Eq. (0) , Eq. (|24|) can be written in the following way, 



2zeukT 



[S]i (zev\ /[S]c / 

[^ ex n2lrJ-\/[S]- exp V 



v/fSUSjlsinh 



ze(v — vg) 



2kT 



Eq. $ 



J V L J^L 1 y ^J, 

is our general result that follows from the combined Ohm's and Fick's law. 



(26) 
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Now the integral / depends upon both the voltage U = U(x) and the cross section A = A(x). To determine the 
x dependence of U(x) one would need Poisson's equation for the electrostatic potential, taking into account the net 
charge distribution in the membrane, including both positive and negative ions. However, this charge distribution 
will depend upon detailed properties of membranes and their channels that have been little known so far. Thus it 
seems a reasonable approach to make certain assumptions directly about U{x). 

A commonly used assumption is that U(x) is linear, i.e. that the electric field —dU/dx is constant, and that the 
cross section is constant, A(x) — Aq. Then Eq. (Eq) takes the form 



Aqv sinh ( - 



ze(v-vs) \ 



As should be expected, this relation simplifies to the usual Ohm's law in the special case [S]i = [S] c , since then v$ = 
by Eq. (Q). Eq. ( p7f ) is known as the Goldman constant field approximation. Goldman (1943) wrote: 

We assume that the membrane contains a large number of dipolar ions near the isotonic point and that these can act to 
minimize distortion in the field especially at low currents. We then approach a situation in which the field is constant and 
are led to a solution analogous to that given by Mott (1939) for electronic conduction in the copper— copper oxide rectifier. 

A more general case, perhaps somewhat more realistic, where the integral / can still be calculated exactly, is that 
of an ion channel having a constant area Aq, except for a short and narrow constriction or pore in its middle, with 
an area A p much smaller than Aq. An example is a cylindrical pore of radius 3 A and length 5 A, which is typical for 
ionic channels. If we furthermore assume a constant electric field everywhere in the channel, and if the length of the 
pore is ed, then we have that 



2dkT \ 1 / zev \ ( 1 1 \ / ezev\ 



Aq \2kTJ \A P A 



The limit of this as v — > is 



In 



Aq A 



ed 

X 



(28) 



(29) 



p 



The last approximation holds when the contribution from the pore dominates the integral, which will be the case 
e.g. when the ratio of areas, A p /Aq, is of the order e 2 . For ev small but nonzero the v dependence of / is only of 
second order in ev, thus it will be a good approximation in a finite voltage range to take / to be constant, equal to Iq. 
The approximation / = constant which is also valid under more general conditions than those assumed in the above 
oversimplified derivation, and it gives 

i = k s sinh ( Z -±^) . (30) 



2kT 

Here kg is independent of v, e.g. in the case considered above 



k s = 2zeukT V[S] e [S]i -± . (31) 

For Na and K ions it is a good approximation to consider the square root of the concentrations ^/[S] e [S]i constant, 
while for Ca ions the relative change in concentration is more significant during one action potential. In the present 
work we used ks = constant in all three cases, for simplicity. We have checked that this does not affect our numerical 
results significantly. 

There is reason to ask whether the linear voltage profile U (x) can be a reasonable approximation in the presence 
of a pore. Indeed, it might seem natural to conclude that most of the voltage drop must be concentrated at the pore 
due to its large resistance. However, with the combined Ohm's and Fick's law, the current is driven by gradients in 
both voltage and concentration, as follows from Eq. (|^). A large current may be due to a large voltage drop over 
the pore, or it may be due to a large change in concentration, without any large voltage difference. Thus, in general 
one has to take into account the detailed properties of the channel in order to see which one of the gradients is the 
dominant driving force in a given situation. 

In a recent investigation by Nonner and Eisenberg (1998), Poisson's equation relating the net charge density and 
electrostatic potential was included in a more extensive analysis for a specific model of a channel with a narrow pore. 
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In their analysis they indeed find that only part of the voltage drop is across the pore (something like half of it). 
In their numerical simulations the voltage in the pore is dominated by the presence of charged carboxyl groups, and 
thus almost independent of the transmembrane voltage. This lends support to the approximation used here, that the 
integral I, Eq. (25), can be regarded as being constant. 

Thus our simple result ( |30| ) has the characteristic features of the current-voltage relationships obtained by Nonner 
and Eisenberg in their more extensive investigation. One characteristic feature is that Eq. (|3(]) shows inward rec- 
tification for large values of [S] /[S]i, i.e. increased conductance for large negative potentials. Indeed the curves in 
figure 3A in Nonner and Eisenberg (1998) have shapes of a hyperbolic sine. Such a behavior is not predicted by the 
Goldman (1943) equation, Eq. (p7|), but is seen in many excitable cells (Hille 1992). This is another reason to base 
our computations on Eq. (|30|) in order to see the consequences of its application. 



3. Potassium Channels 



If the flux of ions is given by Eq. ( |30|) and regulated by the fraction of open channels x, the membrane current 
through potassium channels is 



?k = x sinh 



dx 1 , 
— = — cosh 
at tk 



e(v - v K ) 
2kT 
2e(v - v x ) 
kT 



1 + tanh 



2e(v - f x ) 
kT 



(32) 



where tk = 1/2 A is the maximum value of the relaxation time, fci< is the conductance parameter of Eq. (pl|), vk is 
given by Eq. (|7|), and the time dependence of x is given by Eq. ([ll]) with Eqs. ( |l7| ) and (l2C| ) for Xoo and r respectively. 
Here we have used the identity 



- (1 + tanh( 
2 



1 + exp(-20) 



(33) 



4- Calcium and Sodium Channels 

The calcium and sodium channels have an inactivation mechanism in addition to the above activation mechanism. 
We can view these mechanisms as two independent Markov processes, or equivalently two independent sliding doors, 
which are both affected by voltage. An ion can only go through if both sliding doors are at least slightly open. Here 
the activation mechanism is very fast, with a time constant of only a few milliseconds, so we use the steady state 
fraction of open channels, Eq. (|T7|), for this. The maximum time constant of inactivation for calcium and sodium 
channels are of the same order of magnitude as the maximum time constant of the activation of the potassium channel 
(typically a few hundred milliseconds), thus 



kca f doo sinh 



d °° ~ 2 

£ = J 

dt 



1 + tanh 



e(v - V Ca) 

kT 
2e(v - v d ) 



■ cosh 



TCa 



kT 
2e(v — Vf) 
kf 



(34) 



1 — tanh 



2e(v — Vf) 
kf 



and, 



G 



«Na 

moo 

dh 
~dt 



La h wirTc, sinh 



1 
2 

] 



1 + tanh 



e(v - fNa) 

2kT 
^ 2e(v - v m ) 



■ cosh 



kT 
2e(v - Uh) 
kT 



(35) 



1 — tanh 



2e(i; - v h ) 
kT 



where fcc a and kjsia are the conductance parameters of the calcium and sodium currents respectively, «c a and «Na 
are given by Eqs. (^J) and (^), and u m are the half-activation potentials, and Vf and Vh are the half-inactivation 
potentials. 

Note that the activation and inactivation mechanisms work in the same way, and differ in two respects only. First, 
the time constants differ experimentally by two orders of magnitude, and second, the gating charge q, Eq. ([f4]), is 
+4e in one case and — 4e in the other case. 



C. Na+,K + Pump 

The Na,K~ATPase is found in the plasma membrane of virtually all animal cells and is responsible for active 
transport of sodium and potassium. Low sodium concentration and high potassium concentration in the cytosol are 
essential for basic cellular functions such as excitability, secondary active transport, and volume regulation. In our 
model, the Na + ,K + pump is the only energy source. We shall assume that the following equation is a complete 
macroscopic description of the pump reaction (Chapman 1978), 

ATP + 3Na+ + 2K+ £ ADP + P io + 3Na+ + 2K+ , (36) 

where ATP, ADP and Pj are adenosine triphosphate, adenosine diphosophate and inorganic phosphate, while a and 
(3 are the rates for the forward and backward reactions. The energy involved in the movement of 3 Na + and 2 K + 
ions against their electrochemical gradients is 

AG Na = -3e(v - i) N a) (37) 
AG K - +Mv - «k) , (38) 

where vk and WN a are given by Eqs. (Q) and (^J). This result is independent of the detailed interaction between ions, 
molecules and the ATPase enzyme. Therefore, the total change in Gibbs free energy is 

AG = AGatp + AG Na + AG K 

= e(v A TP + 3u N a - 2v K - v) , (39) 

where AGatp is the energy associated with the breakdown of ATP, and iatp = AGatp/c Note that AG has to be 
negative, at least when averaged over time, but the sum AGNa + AGk may very well be positive, since AGatp is 
large and negative. Thus, part of the energy from ATP breakdown goes into increasing the free energy of potassium 
and sodium ions, but much energy is dissipated, since the energy available is actually much larger than the energy 
required to translocate the potassium and sodium ions at small negative membrane potentials. 

In practice, such a pump or motorized swing door will quickly reach saturation. We therefore choose the sum of the 
forward and backward rates to be constant, resembling the maximum possible speed of the swing door in the forward 
and backward directions, 

a + (3 = X , (40) 

where A is a constant. At equilibrium, the forward reaction must occur just as frequently as the reverse reaction, 
giving 



a 



AG 



p (41) 



Solving Eqs. ( f40| ) and Q) for a and (3 gives 
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The difference 



Aexp(-^) , x 



a - /3 = X f^tMz^m . Atanh (_">) , (44) 
«p(-U)+exp(|f) V *<Tj 



gives the net pump current for a cell with M pumps as 



«NaK = Me(a - (3) = k NaK tanh ( — ) > ( 45 ) 



where &NaK = MeX. 



D. Na + ,Ca 2+ Exchanger 

To maintain a steady state for the intracellular calcium concentration in for example heart cells, the amount of 
calcium that enters the cell via ionic channels must be extruded. The Na + , Ca 2+ exchanger is the major mechanism 
responsible for achieving a balance between calcium entry and extrusion in oscillating cells. We assume that the rate 
for the forward (a) and the backward (j3) exchange reaction given by (Mullins 1977) 

3Na+ + Caf+ A 3Na+ + Ca 2 + , (46) 

are governed largely by the electrochemical gradients for sodium and calcium, together with the membrane potential. 
In other words, the energy produced when 3 extracellular sodium ions take the elevator down into the cytosol is used 
to elevate one calcium ion up from the cytosol into the extracellular space, 

AG Na = +3e(u - v Na ) (47) 
AG Ca = -2e(v - v Ca ) , (48) 

where vc& and UNa are given by Eqs. (^) and (^]). The total work done in reaction ( ff6| ) is 

AG = AG Na + AGca = e(v - 3v Na + 2v Ca ) . (49) 

The ratio of a to j3 in Eq. ( f46|) is again determined by AG like in Eq. (ffl|). However, in the present case saturation 
effects are not expected and furthermore AG will vary around zero, so we put 

/ e(v - 3v Na + 2u Ca )\ ,_„,. 

a=Acx H 2kf — J (50) 

/? = ACX H + 2~kf J' (51) 

where we make the assumption that A is a constant (Mullins, 1981). For a cell with N exchangers the net current is 
then 

«NaCa = -Ne(a - p) = fc NaCa smh I — — I , (52) 

where fc Na Ca = 27VeA. 
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E. Membrane Voltage 



Imagine that the electrical activity of a cell is described by the five currents discussed above, and that all the 

other currents (Boyett 1996) are of minor importance. The standard differential equations for the voltage, and the 
conservation laws for intracellular ionic concentrations are then 

dv 1 . . . . . 

= --jrj (IK + «Ca + «Na + «NaCa + *NaK j , (53) 

|[K]^^f^, (54) 

^ r/-i i 2i NaCa — «Ca /r: r\ 

— [Ca]i = 2FV ' (55) 

|[Na] i= ~ tNa ~ 3 ^~ 3 ' NaCa , (56) 

where C is cell capacitance, F is Faraday's constant, and we assume that the cell volume V is constant. Now Eqs. (|5 
([55|) and (|5^ ) can be solved for ix 7 «Ca 7 and i^ a , and we obtain 

iK = -FV^-[K} i + 2ix aK , (57) 
at 

Z Ca = -2 JV — [Ca]! + 2i N aCa , (58) 

d 

«Na = --FV— [Naji - 3«NaK ~ 3?NaCa ■ (59) 



Inserting this into Eq. (53) yields 



§ = ^|([K]i + 2[Ca] i + [Na] i ) , (60) 



since the remaining currents cancel. This equation can also be written as 



j t (v - ™ {[K], + 2[Ca], + [Na]i}) = . (61) 

This integrated gives 

" - ^ (Mi + 2[Ca]i + [Na]i) = « , (62) 

where the integration constant Vq has to be determined. Given that the voltage across a capacitor is zero when the 
net charge difference is zero, we must choose the integration constant 

FV 

v = -— QK] e + 2[Ca] + [Na] e ) , (63) 

which gives 

FV 

v = — {[K]i - [K] e + 2([Ca]i - [Ca] c ) + [Na]i - [Na] c } . (64) 

This choice of vo depends on the assumption that all other ions have the same concentrations on both sides, consistent 
with Eq. (^||) where it is assumed that they do not contribute to the current. This is also consistent with standard 
assumptions in the literature (Encyclopaedia Britannica 1997) 

In the extracellular fluid, electroncutrality is preserved by a balance between a high concentration of Na+ on the one hand 
and a high concentration of Cl — as well as small quantities of impermeant anions such as bicarbonate, phosphate, and sulfate 
on the other. In the cytoplasm, where K+ concentration is high, the concentration of Cl — is much below that necessary to 
balance the sum of the positive charges. Electroncutrality is maintained there by negatively charged impermeant proteins 
and phosphates. Osmotic balance is maintained between the extracellular fluid and the cytoplasm by movement of water 
through the plasma membrane when the total concentration of particles on one side is not equal to that on the other. 
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Eq. ( p34[) is nothing but the relation between electric potential and charge of a capacitor, which is actually the origin 
of Eq. (p3h. Thus it is completely general and independent of the number of membrane currents in a model. It means 
that: 

The voltage across the membrane of a cell is caused by, and is directly proportional to, the surplus of charge 
inside the cell. 



Since Eq. ( |64| ) is the explicit integral of Eq. (J53| ), it can be used instead of Eq. (53) (or the equivalent of Eq. (p3[)) in 
any model. The differential equation, Eq. fl53[), is needed only in models where the intracellular ionic concentrations 
are not tracked individually (for example, the Hodgkin-Huxley equations (1952)). 

There is a significant difference between Eqs. (|5^) and (|64|) for use in numerical simulations, for the following reason. 
There are two different ways to determine how many ions there are inside a cell. The first method counts every ion 
entering or leaving (Eq. (pf)), while the second method counts all the ions inside the cell (Eq. (|64|)). Both methods 
will give correctly the variation in the number of ions in the cell. However, the observer of ions entering and leaving 
observes only the variations in the number, and if he wants to know the actual number, he must make an initial guess 
of the number of ions already inside. Because his guess may differ significantly from the actual ion number, the results 
from the two methods may be contradictory. 

A variant of Eq. (^) has recently been derived by Varghese and Sell (1997). However, they did not identify the 
integration constant vq, which is related to the initial ionic concentrations and represents the initial net charge via 
the electric capacitance of the cell as shown in Eq. (|64|). 

There is reason to ask whether it is a reasonable approximation to omit the anions in Eq. ([m]). This can be justified 
if the total concentration of cations is approximately the same on both sides. Indeed, this property is seen in most 
ionic models, like for instance in Wilders (1993), where the cation concentrations are 



[K] e = 5.4mM [Ca] c = 2mM [Na] = 140 mM 

[K] ; = 140 mM [Caji = 0.0000804 mM [Naji = 7.5 mM . 



(65) 



F. Energy Balance and Osmotic Pressure 

The current i in Eq. (|2l"l) may be written as 

where n is the number of ions transferred from the inside to the outside of the membrane, and dn/dt is the rate of 
transfer. The change in free energy when one ion is transferred, is AG = — ze(v — vs), and the total change in free 
energy over a time interval is 

pn c-t 

AG = — I ze(v — vs) dn = — I i(v — vs) dt . (67) 
Jo Jo 

It follows from Eq. © that the current i has the same sign as the voltage v — v$ as required in general to have 
thermodynamic stability, so that the integrand in Eq. ( |S7|) is positive (strictly speaking non- negative), and therefore 
AG < (or AG < 0). 

By similar reasoning, taking into account all the reversal potentials and the free energy associated with the break- 
down of ATP, we find that the total change in free energy due to the five currents in our model is 

AG = - / [i K (v - v K ) + i Ca (v - u Ca ) + i Na (u - u Na ) 
Jo 

+ iNaCa(« - 3«Na + 2v C a) + iNaK^ + 2w K - 3«Na ~ VATp)] dt . (68) 

Each of the five terms in the integrand is positive (non- negative) , since each current has the same sign as the 
corresponding voltage. In other words, energy is dissipated all the time by all the five currents, implying that 
AG<0. 

The main contribution to the negative AG is the ATP term, 
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AGATP = / «NaK VATP dt , (69) 

Jo 

which in practice is negative all the time, and furthermore is large in magnitude compared to the other terms. It 
should be noted that this term is the source of useful energy that is dissipated to maintain the activity of the cell 
and keep it away from equilibrium. Keeping this term apart, we may calculate the change in free energy of the ionic 
system, 

AGions = AG — AGatp 

rt 

[i K (v - Vk) + ica(f - «Ca) + «Na(w - l>Na) 



+ «NaCa(« - 3«Na + 2v C a) + ^NaK^ + 2^ K - 3u Na )] dt . (70) 

Using Eqs. (|3|), (|57|), (pq), and ( |59"| ) to eliminate the currents, and assuming the capacitance G and the volume V to 
be constant, we get that 

[K]i 



AG ions - G / vdv-FV v K rf([K]i) 

JO J\K]c 



/•[Cajj ,[Na], 

-2FV vctd([Ca\i) - FV / ^d([Na]i), (71) 

J[Ca] B J[Na] 

where the reversal potentials Vk, ^c a , and vjq a depend on the integration variables [K];, [Ca]j, and [Na]i according to 
Eqs. (0), (^) and Integrating from the equilibrium state v — 0, [K] ; = [K] c , [Ca] ; = [Ca] c , and [Na]i = [Na] e , using 
the indefinite integral 



In 4>d(f) = <t> In 4> -<j>, (72) 



we find 



AG ions = ^ Cv 2 + RTV ([K]j In ( gji )+ [Ca] i In ( + [Na] i In ( [N&]i 



2 T J V[K] C ; L J V[Ca] c ; L J V[Na] c yj (73) 

- RTV ([K]; - [K] c + [Na]; - [Na] c + [Ca] ; - [Ca] c ) . 

In passing it may be noted that in the present case the equilibrium state is simply the one with equal concentrations 
of cations on both sides, as follows from our assumption of having the same concentration of anions or negative charge 
on both sides of the membrane. More generally equilibria for ionic systems are described by the Donnan (1911) 
equilibrium that can yield different concentrations on both sides of the membrane. 

Since AGions is a function only of the state of the cell and is independent of the process by which the state is 
reached, it represents a potential energy for the cell, which we will call P. Note that P = in the equilibrium state, 
whereas P > in all other states. P is the minimum energy needed to bring a thermal system away from equilibrium 
with its surroundings, or equivalently, the maximum work that can be performed by the system when returning to 
equilibrium. 

The potential energy P, as defined in Eq. ([73]), contains three terms, each of which can be given a more direct 
physical interpretation. The first term is simply the electrostatic energy of a capacitor, while the two temperature 
dependent terms are related to thermal properties. In fact, since we assume ideal dilute solutions, the change in 
entropy due to changes of ion concentrations away from their equilibrium values is 

. . UV {[K], 1„ (&) + ,Ca], m (H) + [Na], In (») } . (74, 

Under the same changes, the change in osmotic pressure inside the cell is equal to the difference in osmotic pressure 
across the membrane, which is, for ideal solutions, 

7T = RT ([K]i - [K] c + [Na]; - [Na] + [Ca]i - [Ca] c ) . (75) 

Note that for fixed volume the anions will not contribute to the difference in osmotic pressure, but they will contribute 
if the volume is changed and the membrane is impermeable to them. In terms of the change in transmembrane voltage, 
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v, the change in entropy, s, and the change in transmembrane osmotic pressure, tt, as compared to the equilibrium 
state, we may write 

P = i Cv 2 - Ts - Vir . (76) 



Equation (75) is the van't Hoff equation (1887) for the osmotic pressure across a solute impermeable barrier 
separating two ideal dilute solutions. In 1887 van't Hoff noticed that the behavior of solutes in dilute solutions 
resembles the behavior of a perfect gas (van't Hoff, 1887), and as quoted by Arrhenius in a memoir edited by Jones 
(1899): 

The pressure which a gas exerts at a given temperature if a definite number of molecules is contained in a definite volume, is 
equal to the osmotic pressure which is produced by most substances under the same conditions, if they are dissolved in any 
given liquid. 

Rewriting Eq. (|73|), using Eqs. © and @, we may summarize the energy balance in the following way, 
* 1 

«NaK(v + 2«K - 3t>Na) dt = — Cf 2 sT — TlV 

[i K (v - V K ) + i Ca (v - V C a) + «Na(« ~ «Na) + «NaCa(w ~ 3u N a + 2u C a)] dt. (77) 

The left hand side of this equation is the useful work performed upon the cell by the Na + ,K + pumps, moving Na + 
and K + ions against their potential gradients. The energy supplied by the pumping of ions produces the following 
effects that either change the potential energy of the cell or cause energy loss by dissipation, 

1. a transmembrane voltage difference, v; 

2. a change in entropy, s; 

3. a transmembrane osmotic pressure difference, tt; and 

4. downhill ionic currents through the exchangers and channels, ix, «Ca, *Na, and ZNaCa- 

In an oscillating cell, as described by the present model, the following two inequalities will hold, over a sufficiently 
long time interval, 

- AGatp = / «NaK «atp dt > - i NaK (v + 2v K - 3v Na ) dt > . (78) 
Jo Jo 

The first inequality is simply the inequality ?NaK( u + 2«k — 3«Na — ^atp) > 0, which follows from Eq. J45]). It means 
that the energy released by breakdown of ATP is larger than the useful work performed by the pumps, as required 
from general principles, in other words, that energy is dissipated by the current iNaK produced by the pumps. The 
second inequality must hold due to Eq. (fn\), where the right hand side consists of three oscillating potential energy 
terms plus four positive terms that describe energy dissipation. This inequality shows that the useful work performed 
by the pumps is positive, as required to maintain the dissipation due to the other currents of a working cell away 
from thermal equilibrium. 

We may remark that the laws of Ohm and Fick, equations (pj) and (^), are consistent with the use of ideal 
solutions and osmotic pressure that assume independent (non-interacting) particles. Arrhenius (1902) wrote about 
the relationship between osmotic pressure and diffusion: 

Besides the electrical, other forces may be active in causing the movement of the ions. Of these the osmotic pressure is the 
most important. On account of this pressure a phenomenon called diffusion (hydrodiffusion) may be observed. 

In the model considered the osmotic pressure n has not been involved in the dynamics. However in an extended model 
with variable volume V it will be more important as it will determine the solute flux through the membrane. 
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III. A MODEL FOR CARDIAC PACEMAKER CELLS 



In the above, a mathematical model of the membrane potential has been derived where Eqs . (|7| ), (0), and @ 



represent the equilibrium potentials, Eqs. ( p2[) , (J3J), and ( p5[ ) the ionic currents, Eqs. ( |52| ) and (45) the exchanger 
and the pump currents, Eqs. (^), ( |55| ) and ( |56| ) the ionic concentrations, Eq. ([34]) the membrane voltage, and finally, 
Eq. ( |76"| ) the osmotic pressure across the cell membrane. The model has 6 time dependent variables x, /, h, [K]i, [Ca]i 
and [Na];, and the equations are summarized in Appendix A. 

A. Ionic Mechanisms in the Cardiac Pacemaker 

Akinori Noma published in 1996 an excellent review of the ionic mechanisms of the cardiac pacemaker potential 
(Noma, 1996). In this short paper, Noma investigated the mechanisms that produce spontaneous activity in sinoatrial 
node cells, and introduced the following overview of the relevant ionic currents. 

Channel gating which drives membrane depolarization during diastole 

• Deactivation of «k («Kr)- 

• Removal of inactivation of ica,L and i s t- 

• Activation of the hyperpolarization-activated current (if). 

• Activation of L-type Ca 2+ current (ica,L)- 

• Activation of T-type Ca 2+ current (ica,T)- 



Background conductance 

*t>,Na : A cation current with reversal potential of about — 20 mV. 

*K,ACh : Spontaneous openings of the K + channels. 

*NaK : Na/K pump current. 

*NaCa : Na/Ca exchange current. 

?k,atp : ATP sensitive K + channels. 

We will not try to determine the relative amplitude of the above current components here, but instead demonstrate 
that only five membrane currents is sufficient to ensure stable intracellular ionic concentrations. These are ik, ica, 
*NaK, *NaCa, and iNa- In our model we assume that ica.T and ik s are of minor importance; i.e. when we talk about 
*Ca we mean ica.L, and when we talk about ik we mean i^r- 



B. Model Parameters 



The various parameters play different roles in the model, and we list them in different tables to distinguish between 
fundamental physical constants (table Q), experimentally observed constants (table [il]), adjustable parameters (table 
III) and initial conditions (tabic [Tv| ) in the model. One parameter not listed is the gating charge q, Eq. (|l4|), which 
is the origin of the factor 2e/kT in Eqs. (|32|), (|34|) and (|35|). This corresponds to a slope factor for the activation 
and inactivation curves of kT/Ae « 6.68 mV at 37°C. The observed slope factors are 7.4 mV for activation of zk 
(Shibasaki 1987), 6.6 mV for activation of ic a (Hagiwara et al. 1988), 6.0 mV for inactivation of «c a (Hagiwara et al. 
1988), 6.0mV for activation of i-^ a (Muramatsu et al. 1996), and, finally, 6.4mV for inactivation of ZNa (Muramatsu 
et al. 1996). Hence, we see that kT/Ae, corresponding to a gating charge of q « ±4e, is a good approximation. 

The half-activation and inactivation potentials in the model (i> x , Vd, V{, v m and Vh) are based on the experiments 
of Shibasaki (1987), Hagiwara et al. (1988) and Muramatsu et al. (1996), and we use a value of vatp that gives a 
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reversal potential for the sodium pump in good agreement with the experiments of Sakai et al. (1996). The maximum 
time constants in these experiments were 203 ms for activation of ik (Shibasaki 1987), 225 ms for inactivation of ic a 
(Hagiwara et al. 1988) and 174 ms for inactivation of iNa (Muramatsu et al. 1996) . In the model, however, we combine 
these and use a maximum time constant of 200 ms for both tk, Tc a and TNa- Finally, we use typical values for cell 
volume, cell capacitance, and extracellular ionic concentrations. 

Two of the differential equations in the model have almost identical but opposite dynamics. If we modify the half- 
inactivation potentials of calcium from Vf — —25.0 mV to Vf — v x = —25.1 mV, it is possible to relate the inactivation 
gating of calcium to the activation gating of potassium by the equation 

x + f=l, (79) 

since the time constants for these two processes are equal. We have thus reduced the number of differential equations 
in the model by one. 

This computational saving will be irrelevant for the computation of one action potential, but is important in an 
extended model with thousands of coupled cell, or in a long-time integration of the one cell model. 



C. Pacemaker Current 



The relative amplitude of the ionic currents that drive membrane depolarization during diastole is still a matter 
of debate. DiFrancesco (1993) argues that the hyperpolarization activated current (if) is the only current that can 
generate and control the slow depolarization of pacemaker cells. The if current is normally carried by Na + and K + . 
Guo et al. (1995a; 1995b; 1996) reported another current, called the sustained inward current i s t, where the major 
charge carrier is believed to be Na + . Also a Ca 2+ "window" current has been observed in rabbit sinoatrial node cells 
(Denyer & Brown 1990). It is possible that any one of these currents, or a combination of them, is responsible for 
membrane depolarization during diastole. However, the estimates of the net membrane current during diastole is so 
imprecise (Zaza et al. 1997), that we could not form a judgment on the question. 

During diastole the electrochemical driving forces produce outward K + currents and inward Na + and Ca 2+ currents, 
and the driving force for Ca 2+ is much larger than that for Na + . These findings implies that a significant background 
influx of Ca 2+ is possible during diastole, and that this current might be responsible for the pacemaker activity in 
sinoatrial node cells. We denote the conductance for this current &b,Ca, an d modify (|1|) combined with ( f79| ) to read 
as 

, (80) 

V T J 

with vt — kT/e. In our model /cb.Ca is responsible for the slow diastolic depolarization, and it is thus our pacemaker 
current. However, as suggested by Guo et al. (1995b), it is indeed possible that the sustained inward current i s t may 
largely replace the role of the Ca 2+ currents, assumed here and in previous studies (Wilders 1993). 



D. Adjustable Parameters 



The density of ionic channels, exchangers and pumps (i.e. fcc a , fcb.Ca, &Na, &k, fosiaK and &NaCa) can vary significantly 
from cell to cell. In order to reproduce recorded action potentials (Fig. 7 A. in Baruscotti et al. (1996)), we fit the 
adjustable parameters (table til) and the initial conditions (table IV) numerically. More details of the method are 
given in (Endresen 1997a). Many different combinations of fco, &b,Ca, ^ Na ' ^NaK, and &NaCa resulted in good 
approximations to the experimentally recorded waveform, from which we conclude that different cells can produce 
the same action potential although they have different mixtures of ionic channels, exchangers, and pumps. 



E. Simulation Results 



The five differential equations in the model were solved numerically using a fifth-order Rungc-Kutta method with 
variable steplength. More details are given in (Endresen 1997b). We computed the work W, defined as minus the 
integral on the right hand side of (f70|), to check that the equation W + P = was satisfied numerically. This could also 
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be used for varying the steplength in an efficient way (Marthinsen et al. 1997), since the solution of our differential 
equations must satisfy this constraint. These "checksum equations" are shown in Appendix B. 

In Fig. 1 (a) the modeled action potential is shown together with the experimental curve of Baruscotti et al. (1996). 
The curves are identical in shape, but we adjusted the modeled curve somewhat (we multiplied the voltage amplitude 
by a factor 1.25, without changing the minimum value) to obtain the same voltage amplitudes. At present it is not 
clear to us which mechanism is needed in the model in order to avoid this factor. Fig. 1 (b) shows the five membrane 
currents in the model, z Ca , «Na, «k, «NaK, and ^NaCa- 

Fig. 2 shows the spontaneous action potentials together with the intracellular ionic concentrations and the osmotic 
pressure across the cell membrane. These computations used the initial conditions stated in table IV. Cells must 
generate their membrane potential by actively transporting ions against the respective concentration gradients. To 
examine this process in our model, we ran a simulation starting with equal intracellular and extracellular ionic 
concentrations: [K] ; = [K] c = 5.4 mM, [Ca]i = [Ca] c = 2mM, and [Na]i = [Na] c = 140 mM. The results are presented 
in Fig. 3, that shows the voltage and Nernst potentials (a), and the energies (b) in a long time simulation. After 
approximately 750 seconds (12.5 minutes) the system reaches oscillations identical to the original oscillations shown 
in Figs. 1 and 2 (this can not be seen from Fig. 3 since the time scale is very different). This long time simulation is 
a numerical indication that the oscillations in Fig. 2 and 3 indeed correspond to a stable limit cycle. 



IV. DISCUSSION 

We have presented a simple model for the cells of the rabbit sinoatrial node. Our model involves only Na + , K + , and 
Ca 2+ ions, their respective channels, the Na + ,Ca 2+ exchanger, and the Na + ,K + pump. The equations were derived 
using basic physical principles and conservation laws. Since the only source of energy in our model is the sodium 
potassium pump, we can easily track the flow of energy, and show that the pump works to generate a transmembrane 
voltage, osmotic pressure difference, and an entropy. Our equations also account for the energy lost due to downhill 
ionic fluxes through the exchanger and channels. A prediction of osmotic pressure variations is a novel result of our 
energy analysis. 

The intracellular ionic concentrations are dynamic variables in our model, governed by the conservation Eqs. (p4), 



fl55p, and (56). This allows us to replace the standard differential equation for the voltage (|53J) with the algebraic 
Eq. (pj3) . Although a number of other ionic models also keep track of intracellular ionic concentrations (see Wilders 
(1993)), we are unaware of any other model using an algebraic equation for the membrane potential. Models that use 
the standard voltage differential Eq. ( |5"3| ) have a phase space with one superfluous extra dimension. The initial condi- 
tion for this extra differential equation cannot be chosen independently of the initial conditions for the conservation 
Eqs. (|54|), (|5"5|), and ( |56| ) - if it is, the computed membrane potential will be erroneous. For these reasons, we suggest 
that our algebraic expression for the membrane potential should replace the standard voltage differential equation in 
models where intracellular ionic concentrations are dynamic variables. 

Our model does not include the funny current (if), ATP sensitive channels, stretch-activated channels, or other ion 
channels that may be important (Boyett 1996). We also ignored the effect of calcium uptake and release from the sar- 
coplasmatic reticulum, which would affect the Nernst potential of calcium, but not the membrane potential. We have 
assumed that the ionic channels arc governed by a Markov process, that the maximum of the activation/inactivation 
time constant occurs at the same voltage as the inflection point of the sigmoidal steady state activation/inactivation 
curve, and that the steady state activation/inactivation curves were temperature independent. Also, we have assumed 
that the cell volume is constant. While such assumptions reduce the number of parameters in the model, they may 
also result in discrepancies with experiment. 

Finally, we would like to point out that our model is based on experiments where some were conducted at room 
temperature (22-24°C) (Baruscotti et al. 1996; Muramatsu et al. 1996), while others were performed at 37°C 
(Shibasaki 1987; Hagiwara et al. 1988; Sakai et al., 1996). It is not clear what effect varying temperature has in our 
model as this was not checked out numerically. 



The values of the parameters &ca, &Na, fe, &NaK and /cNaCa, given in table III, arc only an estimate of the actual 
physiological parameters. We did not systematically study the dynamics of the model for different values of the 
parameters, but we hope that future experiments will help to discriminate between different parameter sets that may 
reproduce the experimentally recorded action potentials. 
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APPENDIX A: EQUATIONS OF MOTION 
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APPENDIX B: CHECKSUM EQUATION: W + P = 

i K (v - V K ) + ica(v - V Ca .) + *Na(« - "Na) 

+»NaCa(w - 3u Na + 2v C a) + ^NaK^ + 2« K - 3« Na ) 

1 
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TABLE I. Fundamental Physical Constants 



Name 


Value 


Unit 


k 


1.38065812 • 10 _2U 


mJ/K 


e 


1.60217733 • 10 -19 


C 


F 


96485.30929 


C/mol 


R = kF/e 


8314.511935 


J/kmol K 


TABLE II. Observed Constants 


Name 


Value 


Unit 


T 


310.15 


K 


[K]e 


5.4 


mM 


[Ca] e 


2 


mM 


[Na] e 


140 


mM 


V 


10 


10 3 /im 3 


C 


47 


pF 


fx 


-25.1 


mV 


Vd 


-6.6 


mV 


Vf 


-25.0 


mV 


Dm 


-41.4 


mV 


Vh 


-91.0 


mV 


VATP 


-450 


mV 


T — TK = TCa = TNa 


200 


ms 


v T = kT/e = RT/F 


26.7268 


mV 


TABLE III. Adjustable Parameters 


Name 


Value 


Unit 


fcca 


26.2 


pA 


feb.Ca 


0.01645 


pA 




112.7 


pA 




32.9 


pA 


feNaCa 


1400.0 


pA 


^NaK 


11.46 


pA 


TABLE IV. Initial Conditions 


Name 


Value 


Unit 




0.1 




fo = l-x 


0.9 




h 


0.008 




[K]i 


130.66 


mM 


[Ca] i0 


0.0006 


mM 


[Na] i0 


18.7362 


mM 
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FIG. 1. Action potential and currents, (a) Experimentally recorded and scaled (by a factor 1.25) model-generated rabbit 
sinoatrial action potential waveform, (b) The outward delayed rectifying potassium current (jk), the inward calcium current 
(ica), the inward sodium current (iNa), the sodium calcium exchange current (iNaCa) and the sodium potassium pump current 
(*NaK). These computations used the initial conditions in table £y| 
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FIG. 2. Membrane potential (not scaled), intracellular ionic concentrations and osmotic pressure of a rabbit sinoatrial node 
cell, (a) Model-generated action potential waveform, (b) potassium concentration [K]i, (c) calcium concentration [Ca]i, (d) 
sodium concentration [Na]j and (e) the osmotic pressure tt across the cell membrane. These computations used the initial 
conditions in table 0. 
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FIG. 3. Long time simulation showing the membrane potential, Nernst potentials and energies starting with equal intracel- 
lular and extracellular concentrations: [K] ; = [K] c = 5.4 mM, [Ca]i = [Ca] c = 2mM and [Na] ; = [Na] e = 140 mM. (a) Nernst 
potential for calcium (wca), potassium («k), sodium (t)N»), and membrane potential v. (b) Work (W), potential energy (P) 
and total energy balance (W + P). 
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